##Goals

Load in phyloseq data with rooted tree. Evaluate sequencing depth and remove sample. Normalize the read counts between samples. Calculate community dissimilarities. Numbers between 0 and 1. If 0, completely similar versus if they are 1, then they’re completely dissimilar. Sorensen: Shared Species as a binary value: Abundance-unweighted Bray-Curtis: Shared Abundant species: Abundance-weighted (Abundance-)Weighted UNIFRAC: Consider Abundant Species and where they fall on the tree Visualize the community data with two unconstrained Ordinations: PCoA: Linear Method. Eigenvalue = how much variation is explained by each axis. Choose to view axis 1, 2, 3, etc. and plot them together. NMDS: Non-linear. Smush multiple Dimensions into 2 or 3 axes. Need to report Stress value (ideally <0.15). Run statistics with PERMANOVA and betadispR.

##Set seed

set.seed(2443)

##Load libraries

#install.packages("vegan")
pacman::p_load(tidyverse, devtools, phyloseq, patchwork, vegan,
               install = FALSE)
library(gridExtra)
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine
# Load Lakestie colors 
lakesite_colors <- c("110" = "dodgerblue4","M15" = "dodgerblue2","M45" = "#D9CC3C","MLB" = "#A0E0BA")

##Load Data

# Load in rooted phylogenetic tree! 
load("data/03_Phylogenetic_Tree/phytree_preprocessed_physeq.RData")

##Explore Read Counts

# Calculate the total number of reads per sample. 
raw_TotalSeqs_df_combined  <- 
  midroot_physeq %>%
  # calculate the sample read sums 
  sample_sums() %>%
  data.frame()
# name the column 
colnames(raw_TotalSeqs_df_combined)[1] <- "TotalSeqs"
  
raw_TotalSeqs_df_DNA  <- 
  midroot_physeq_DNA %>%
  # calculate the sample read sums 
  sample_sums() %>%
  data.frame()
# name the column 
colnames(raw_TotalSeqs_df_DNA)[1] <- "TotalSeqs"
  
raw_TotalSeqs_df_RNA  <- 
  midroot_physeq_RNA %>%
  # calculate the sample read sums 
  sample_sums() %>%
  data.frame()
# name the column 
colnames(raw_TotalSeqs_df_RNA)[1] <- "TotalSeqs"
  

head(raw_TotalSeqs_df_combined)
##            TotalSeqs
## 110D2W115D     13235
## 110D2W115R     12953
## 110D2W415D     14765
## 110D2W415R     11587
## 110D2W515D     18867
## 110D2W515R     12432
head(raw_TotalSeqs_df_DNA)
##            TotalSeqs
## 110D2W115D     13235
## 110D2W415D     14765
## 110D2W515D     18867
## 110D2W615D      9289
## 110D2W715D     14757
## 110D2W815D     36720
head(raw_TotalSeqs_df_RNA)
##            TotalSeqs
## 110D2W115R     12953
## 110D2W415R     11587
## 110D2W515R     12432
## 110D2W615R     16887
## 110D2W715R     14711
## 110D2W815R     10574
# make histogram of raw reads 
raw_TotalSeqs_df_combined %>%
  ggplot(aes(x = TotalSeqs)) + 
  geom_histogram(bins = 50) + 
  scale_x_continuous(limits = c(0, 10000)) + 
  labs(title = "Raw Combined Sequencing Depth Distribution") + 
  theme_classic()
## Warning: Removed 74 rows containing non-finite outside the scale range
## (`stat_bin()`).
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_bar()`).

raw_TotalSeqs_df_DNA %>%
  ggplot(aes(x = TotalSeqs)) + 
  geom_histogram(bins = 50) + 
  scale_x_continuous(limits = c(0, 10000)) + 
  labs(title = "Raw DNA Sequencing Depth Distribution") + 
  theme_classic()
## Warning: Removed 43 rows containing non-finite outside the scale range (`stat_bin()`).
## Removed 2 rows containing missing values or values outside the scale range
## (`geom_bar()`).

raw_TotalSeqs_df_RNA %>%
  ggplot(aes(x = TotalSeqs)) + 
  geom_histogram(bins = 50) + 
  scale_x_continuous(limits = c(0, 10000)) + 
  labs(title = "Raw RNA Sequencing Depth Distribution") + 
  theme_classic()
## Warning: Removed 31 rows containing non-finite outside the scale range (`stat_bin()`).
## Removed 2 rows containing missing values or values outside the scale range
## (`geom_bar()`).

#What are the minimum number of sequences for each rooted physeq object?
midroot_physeq %>% sample_sums() %>% min()
## [1] 2142
midroot_physeq_DNA %>% sample_sums() %>% min()
## [1] 2582
midroot_physeq_RNA %>% sample_sums() %>% min()
## [1] 2142
#The minimum number of sequences for DNA samples = 2582 and for RNA samples = 2142
#In order to normailze the read counts we first need to import a scale read function from http://deneflab.github.io/MicrobeMiseq/

matround <- function(x){trunc(x+0.5)}

scale_reads <- function(physeq, n = min(sample_sums(physeq)), round = "round") {
  
  # transform counts to n
  physeq.scale <- transform_sample_counts(physeq, function(x) {(n * x/sum(x))})
  
  # Pick the rounding functions
  if (round == "floor"){
    otu_table(physeq.scale) <- floor(otu_table(physeq.scale))
  } else if (round == "round"){
    otu_table(physeq.scale) <- round(otu_table(physeq.scale))
  } else if (round == "matround"){
    otu_table(physeq.scale) <- matround(otu_table(physeq.scale))
  }
  
  # Prune taxa and return new phyloseq object
  physeq.scale <- prune_taxa(taxa_sums(physeq.scale) > 0, physeq.scale)
  return(physeq.scale)
}
#Scale the reads using the above function, not we are scaling each object by their minimum number for sequences.

# Scale reads combined DNA/RNA by the above function
scaled_rooted_physeq <- 
   midroot_physeq %>%
  scale_reads(round = "matround")

# Scale reads DNA samples by the above function
scaled_rooted_physeq_DNA <- 
  midroot_physeq_DNA %>%
  scale_reads(round = "matround")

# Scale reads RNA samples by the above function
scaled_rooted_physeq_RNA <- 
  midroot_physeq_RNA %>%
  scale_reads(round = "matround")

# Calculate combined DNA//RNA read depth 
scaled_TotalSeqs_df <- 
  scaled_rooted_physeq %>%
  sample_sums() %>%
  data.frame()
colnames(scaled_TotalSeqs_df)[1] <-"TotalSeqs"

head(scaled_TotalSeqs_df) #Inspect to see if the scaling occured
##            TotalSeqs
## 110D2W115D      2143
## 110D2W115R      2130
## 110D2W415D      2137
## 110D2W415R      2140
## 110D2W515D      2146
## 110D2W515R      2147
# Calculate DNA samples read depth 
scaled_TotalSeqs_DNA_df <- 
  scaled_rooted_physeq_DNA %>%
  sample_sums() %>%
  data.frame()
colnames(scaled_TotalSeqs_DNA_df)[1] <-"TotalSeqs"

head(scaled_TotalSeqs_DNA_df) #Inspect to see if the scaling occured
##            TotalSeqs
## 110D2W115D      2581
## 110D2W415D      2583
## 110D2W515D      2584
## 110D2W615D      2586
## 110D2W715D      2583
## 110D2W815D      2568
# Calculate RNA samples read depth 
scaled_TotalSeqs_RNA_df <- 
  scaled_rooted_physeq_RNA %>%
  sample_sums() %>%
  data.frame()
colnames(scaled_TotalSeqs_RNA_df)[1] <-"TotalSeqs"

# Inspect
head(scaled_TotalSeqs_RNA_df)
##            TotalSeqs
## 110D2W115R      2130
## 110D2W415R      2140
## 110D2W515R      2147
## 110D2W615R      2138
## 110D2W715R      2142
## 110D2W815R      2135
#Check the range of the data 
#For Combined,DNA and RNA
min_seqs_combined <- min(scaled_TotalSeqs_df$TotalSeqs);min_seqs_combined
## [1] 2117
min_seqs_DNA <- min(scaled_TotalSeqs_DNA_df$TotalSeqs); min_seqs_DNA
## [1] 2568
min_seqs_RNA <- min(scaled_TotalSeqs_RNA_df$TotalSeqs); min_seqs_RNA
## [1] 2129
max_seqs_combined <- max(scaled_TotalSeqs_df$TotalSeqs);max_seqs_combined
## [1] 2173
max_seqs_DNA <- max(scaled_TotalSeqs_DNA_df$TotalSeqs); max_seqs_DNA
## [1] 2595
max_seqs_RNA <- max(scaled_TotalSeqs_RNA_df$TotalSeqs); max_seqs_RNA
## [1] 2173
#Range 
max_seqs_combined - min_seqs_combined
## [1] 56
max_seqs_DNA - min_seqs_DNA
## [1] 27
max_seqs_RNA - min_seqs_RNA
## [1] 44
# Plot Histogram 
scaled_TotalSeqs_df %>%
  ggplot(aes(x = TotalSeqs)) + 
  geom_histogram(bins = 50) + 
  scale_x_continuous(limits = c(0, 10000)) + 
  labs(title = "Scaled Sequencing Depth at 2142") + 
  theme_classic()
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_bar()`).

scaled_TotalSeqs_DNA_df %>%
  ggplot(aes(x = TotalSeqs)) + 
  geom_histogram(bins = 50) + 
  scale_x_continuous(limits = c(0, 10000)) + 
  labs(title = "Scaled Sequencing Depth of DNA only samples at 2582") + 
  theme_classic()
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_bar()`).

scaled_TotalSeqs_RNA_df %>%
  ggplot(aes(x = TotalSeqs)) + 
  geom_histogram(bins = 50) + 
  scale_x_continuous(limits = c(0, 10000)) + 
  labs(title = "Scaled Sequencing Depth of RNA only samples at 2142") + 
  theme_classic()
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_bar()`).

#We can now start running PCoA with respect to Soreson and Bray-Curtis

# Calculate sorensen dissimilarity: Abundance-unweighted of shared taxa
# Running PCoA with limnion as the shape 
scaled_soren_pcoa_DNA <-  
  ordinate(
    physeq = scaled_rooted_physeq_DNA,
    method = "PCoA",
    distance = "bray", binary = TRUE)

#str(scaled_soren_pcoa_DNA)

# Plot the ordination 
soren_limnion_pcoa_DNA <- plot_ordination(
  physeq = scaled_rooted_physeq_DNA,
  ordination = scaled_soren_pcoa_DNA,
  color = "lakesite",
  shape = "limnion",
  title = "Sorensen PCoA All Taxa (DNA)") +
  geom_point(size=5, alpha = 0.5, aes(color = lakesite)) +
  scale_shape_manual(values = c(15,16,17,18, 19, 20)) +
  scale_color_manual(values = lakesite_colors) + 
  theme_bw()

# Calculate sorensen dissimilarity: Abundance-unweighted of shared taxa
scaled_soren_pcoa_RNA <-  
  ordinate(
    physeq = scaled_rooted_physeq_RNA,
    method = "PCoA",
    distance = "bray", binary = TRUE)

#str(scaled_soren_pcoa_DNA)

# Plot the ordination 
soren_limnion_pcoa_RNA <- plot_ordination(
  physeq = scaled_rooted_physeq_RNA,
  ordination = scaled_soren_pcoa_RNA,
  color = "lakesite",
  shape = "limnion",
  title = "Sorensen PCoA Active Taxa (RNA)") +
  geom_point(size=5, alpha = 0.5, aes(color = lakesite)) +
  scale_shape_manual(values = c(15,16,17,18, 19, 20)) +
  scale_color_manual(values = lakesite_colors) + 
  theme_bw()

# Show the plot 
grid.arrange(soren_limnion_pcoa_DNA,soren_limnion_pcoa_RNA, ncol=2) 

#Running PCoA with season as the shape 
# Plot the ordination 
soren_season_pcoa_DNA <- plot_ordination(
  physeq = scaled_rooted_physeq_DNA,
  ordination = scaled_soren_pcoa_DNA,
  color = "lakesite",
  shape = "season",
  title = "Sorensen PCoA DNA") +
  geom_point(size=5, alpha = 0.5, aes(color = lakesite)) +
  scale_shape_manual(values = c(15,16,17,18, 19, 20)) +
  scale_color_manual(values = lakesite_colors) + 
  theme_bw()



# Plot the ordination 
soren_season_pcoa_RNA <- plot_ordination(
  physeq = scaled_rooted_physeq_RNA,
  ordination = scaled_soren_pcoa_RNA,
  color = "lakesite",
  shape = "season",
  title = "Sorensen PCoA RNA") +
  geom_point(size=5, alpha = 0.5, aes(color = lakesite)) +
  scale_shape_manual(values = c(15,16,17,18, 19, 20)) +
  scale_color_manual(values = lakesite_colors) + 
  theme_bw()

# Show the plot 
grid.arrange(soren_season_pcoa_DNA,soren_season_pcoa_RNA,ncol=2)

The PCoAs explain a fair amount of variation for RNA/DNA about 47% for DNA and 40.3% for RNA. The clustering seems quite similar between the DNA and RNA samples this indicates with the MLB(surface) samples seem to cluster in two regions. The liminions do not seem to match indicating that the the variation may not be due to the liminion. The other samples seem to be spread quite wide with no clear pattern arising. Though the DNA and RNA seem to cluster similarly a point of note is that there seems to be a larger difference in the RNA surface MLB cluster on the second axis indicating that there is a difference in active taxa relative to total taxa, this does not seem to be limnion dependant.

# Calculate the BC distance
scaled_BC_pcoa_DNA <- 
  ordinate(
    physeq = scaled_rooted_physeq_DNA,
    method = "PCoA",
    distance = "bray")

scaled_BC_pcoa_RNA <- 
  ordinate(
    physeq = scaled_rooted_physeq_RNA,
    method = "PCoA",
    distance = "bray")



# Plot the PCoA
bray_pcoa_limnion_DNA <- 
  plot_ordination(
    physeq = scaled_rooted_physeq_DNA,
    ordination = scaled_BC_pcoa_DNA,
    color = "lakesite",
    shape = "limnion",
    title = "Bray-Curtis PCoA DNA Limnion") +
  geom_point(size=5, alpha = 0.5, aes(color = lakesite)) +
  scale_shape_manual(values = c(15,16,17,18, 19, 20,21)) +
  scale_color_manual(values = c(lakesite_colors)) + 
  theme_bw()

bray_pcoa_limnion_RNA <- 
  plot_ordination(
    physeq = scaled_rooted_physeq_RNA,
    ordination = scaled_BC_pcoa_RNA,
    color = "lakesite",
    shape = "limnion",
    title = "Bray-Curtis PCoA RNA Limnion") +
  geom_point(size=5, alpha = 0.5, aes(color = lakesite)) +
  scale_shape_manual(values = c(15,16,17,18, 19, 20,21)) +
  scale_color_manual(values = c(lakesite_colors)) + 
  theme_bw()

bray_pcoa_season_DNA <- 
  plot_ordination(
    physeq = scaled_rooted_physeq_DNA,
    ordination = scaled_BC_pcoa_DNA,
    color = "lakesite",
    shape= "season",
    title = "Bray-Curtis PCoA DNA Season") +
  geom_point(size=5, alpha = 0.5, aes(color = lakesite)) +
  scale_shape_manual(values = c(15,16,17,18, 19, 20,21)) +
  scale_color_manual(values = c(lakesite_colors)) + 
  theme_bw()

bray_pcoa_season_RNA <- 
  plot_ordination(
    physeq = scaled_rooted_physeq_RNA,
    ordination = scaled_BC_pcoa_RNA,
    color = "lakesite",
    shape= "season",
    title = "Bray-Curtis PCoA RNA Season") +
  geom_point(size=5, alpha = 0.5, aes(color = lakesite)) +
  scale_shape_manual(values = c(15,16,17,18, 19, 20,21)) +
  scale_color_manual(values = c(lakesite_colors)) + 
  theme_bw()



#Show limnion plot
grid.arrange(bray_pcoa_limnion_DNA,
bray_pcoa_limnion_RNA, ncol=2)

#Show season plot
grid.arrange(bray_pcoa_season_DNA,
bray_pcoa_season_RNA, ncol=2) 

The PCOAs explain a high percentage of the variance for both the DNA (~45%) and RNA(~40%) samples. The top plots look at the limnion dependence of the abundant taxa and we can see that the points aren’t clustering based on shape. However, we see that the points cluster closely based on both shape (representing season) and color (representing lakesite) for the season plots. We see that the surface taxa once again cluster distinctly along the first axis for all these plots, indicating that they are quite different from the other depths that were sampled.

Weighted Unifrac

scaled_wUNI_pcoa_combined <- 
  ordinate(
    physeq = scaled_rooted_physeq,
    method = "PCoA",
    distance = "wunifrac")

scaled_wUNI_pcoa_DNA <- 
  ordinate(
    physeq = scaled_rooted_physeq_DNA,
    method = "PCoA",
    distance = "wunifrac")

scaled_wUNI_pcoa_RNA <- 
  ordinate(
    physeq = scaled_rooted_physeq_RNA,
    method = "PCoA",
    distance = "wunifrac")


# Plot the PCoA
wUNI_station_pcoa_combined_season <- 
  plot_ordination(
    physeq = scaled_rooted_physeq,
    ordination = scaled_wUNI_pcoa_combined,
    color = "lakesite",
    shape = "season",
    title = "Weighted Unifrac PCoA Combined season") +
  geom_point(size=5, alpha = 0.5, aes(color = lakesite)) +
  scale_shape_manual(values = c(15,16,17,18, 19, 20)) +
  scale_color_manual(values = c(lakesite_colors)) + 
  theme_bw()

wUNI_station_pcoa_combined_limnion <- 
  plot_ordination(
    physeq = scaled_rooted_physeq,
    ordination = scaled_wUNI_pcoa_combined,
    color = "lakesite",
    shape = "limnion",
    title = "Weighted Unifrac PCoA Combined limnion") +
  geom_point(size=5, alpha = 0.5, aes(color = lakesite)) +
  scale_shape_manual(values = c(15,16,17,18, 19, 20)) +
  scale_color_manual(values = c(lakesite_colors)) + 
  theme_bw()


wUNI_station_pcoa_DNA_season <- 
  plot_ordination(
    physeq = scaled_rooted_physeq_DNA,
    ordination = scaled_wUNI_pcoa_DNA,
    color = "lakesite",
    shape = "season",
    title = "Weighted Unifrac PCoA DNA season") +
  geom_point(size=5, alpha = 0.5, aes(color = lakesite)) +
  scale_shape_manual(values = c(15,16,17,18, 19, 20)) +
  scale_color_manual(values = c(lakesite_colors)) + 
  theme_bw()

wUNI_station_pcoa_DNA_limnion <- 
  plot_ordination(
    physeq = scaled_rooted_physeq_DNA,
    ordination = scaled_wUNI_pcoa_DNA,
    color = "lakesite",
    shape = "limnion",
    title = "Weighted Unifrac PCoA DNA limnion") +
  geom_point(size=5, alpha = 0.5, aes(color = lakesite)) +
  scale_shape_manual(values = c(15,16,17,18, 19, 20)) +
  scale_color_manual(values = c(lakesite_colors)) + 
  theme_bw()

wUNI_station_pcoa_RNA_season <- 
  plot_ordination(
    physeq = scaled_rooted_physeq_RNA,
    ordination = scaled_wUNI_pcoa_RNA,
    color = "lakesite",
    shape = "season",
    title = "Weighted Unifrac PCoA RNA season") +
  geom_point(size=5, alpha = 0.5, aes(color = lakesite)) +
  scale_shape_manual(values = c(15,16,17,18, 19, 20)) +
  scale_color_manual(values = c(lakesite_colors)) + 
  theme_bw()

wUNI_station_pcoa_RNA_limnion <- 
  plot_ordination(
    physeq = scaled_rooted_physeq_RNA,
    ordination = scaled_wUNI_pcoa_RNA,
    color = "lakesite",
    shape = "limnion",
    title = "Weighted Unifrac PCoA RNA limnion") +
  geom_point(size=5, alpha = 0.5, aes(color = lakesite)) +
  scale_shape_manual(values = c(15,16,17,18, 19, 20)) +
  scale_color_manual(values = c(lakesite_colors)) + 
  theme_bw()

grid.arrange(wUNI_station_pcoa_combined_season,wUNI_station_pcoa_combined_limnion,ncol=2)

#Difference between seasons
grid.arrange(wUNI_station_pcoa_DNA_season, wUNI_station_pcoa_RNA_season,ncol=2)

#Difference between depths
grid.arrange(wUNI_station_pcoa_DNA_limnion, wUNI_station_pcoa_RNA_limnion,ncol=2)

The PCOAs explain a high percentage of the variance for both the DNA (~45%) and RNA(~47%) samples. We do not see a limnion or season dependence for the weighted uniFRAC distance. We did note one interesting feature, however. For the 110 site, there seems to be a cluster in the DNA samples showing that the total taxa present at this site are quite evolutionary similar. The active taxa, however, which we are inferring from the RNA samples show no clustering indicating that they are phylogenetically different.

# Calculate the Weighted Unifrac distance
scaled_wUNI_nmds_combined <- 
  ordinate(
    physeq = scaled_rooted_physeq,
    method = "NMDS",
    distance = "wunifrac")
## Run 0 stress 0.1418455 
## Run 1 stress 0.14491 
## Run 2 stress 0.1615657 
## Run 3 stress 0.1562387 
## Run 4 stress 0.166838 
## Run 5 stress 0.1619733 
## Run 6 stress 0.1494274 
## Run 7 stress 0.1477715 
## Run 8 stress 0.1550695 
## Run 9 stress 0.147955 
## Run 10 stress 0.1542195 
## Run 11 stress 0.1449249 
## Run 12 stress 0.1542195 
## Run 13 stress 0.1668388 
## Run 14 stress 0.1449098 
## Run 15 stress 0.1418455 
## ... New best solution
## ... Procrustes: rmse 1.651949e-05  max resid 0.0001469256 
## ... Similar to previous best
## Run 16 stress 0.1634479 
## Run 17 stress 0.1604705 
## Run 18 stress 0.147952 
## Run 19 stress 0.1600244 
## Run 20 stress 0.1418454 
## ... New best solution
## ... Procrustes: rmse 3.563792e-05  max resid 0.0003311476 
## ... Similar to previous best
## *** Best solution repeated 1 times
scaled_wUNI_nmds_DNA <- 
  ordinate(
    physeq = scaled_rooted_physeq_DNA,
    method = "NMDS",
    distance = "wunifrac")
## Run 0 stress 0.1925474 
## Run 1 stress 0.1985646 
## Run 2 stress 0.2141103 
## Run 3 stress 0.1985623 
## Run 4 stress 0.1925474 
## ... Procrustes: rmse 4.622689e-06  max resid 1.70631e-05 
## ... Similar to previous best
## Run 5 stress 0.2515052 
## Run 6 stress 0.1925474 
## ... Procrustes: rmse 1.281947e-05  max resid 8.138397e-05 
## ... Similar to previous best
## Run 7 stress 0.1926273 
## ... Procrustes: rmse 0.004161513  max resid 0.02407351 
## Run 8 stress 0.1925474 
## ... Procrustes: rmse 3.626435e-05  max resid 0.0002383358 
## ... Similar to previous best
## Run 9 stress 0.2102099 
## Run 10 stress 0.1926274 
## ... Procrustes: rmse 0.004159529  max resid 0.02406611 
## Run 11 stress 0.198565 
## Run 12 stress 0.1985652 
## Run 13 stress 0.2101843 
## Run 14 stress 0.1926273 
## ... Procrustes: rmse 0.004166649  max resid 0.02412437 
## Run 15 stress 0.1925474 
## ... Procrustes: rmse 1.676562e-05  max resid 0.0001045908 
## ... Similar to previous best
## Run 16 stress 0.1985644 
## Run 17 stress 0.1985654 
## Run 18 stress 0.1925474 
## ... Procrustes: rmse 6.392559e-06  max resid 2.872485e-05 
## ... Similar to previous best
## Run 19 stress 0.2451172 
## Run 20 stress 0.1925474 
## ... Procrustes: rmse 7.540351e-06  max resid 4.720985e-05 
## ... Similar to previous best
## *** Best solution repeated 6 times
scaled_wUNI_nmds_RNA <- 
  ordinate(
    physeq = scaled_rooted_physeq_RNA,
    method = "NMDS",
    distance = "wunifrac")
## Run 0 stress 0.1645868 
## Run 1 stress 0.2660364 
## Run 2 stress 0.1645868 
## ... New best solution
## ... Procrustes: rmse 4.067466e-06  max resid 2.252153e-05 
## ... Similar to previous best
## Run 3 stress 0.1645868 
## ... Procrustes: rmse 4.083904e-06  max resid 1.905353e-05 
## ... Similar to previous best
## Run 4 stress 0.1645868 
## ... Procrustes: rmse 2.722264e-06  max resid 1.323778e-05 
## ... Similar to previous best
## Run 5 stress 0.1645868 
## ... Procrustes: rmse 1.98128e-06  max resid 1.034146e-05 
## ... Similar to previous best
## Run 6 stress 0.1645868 
## ... Procrustes: rmse 2.849117e-06  max resid 1.05118e-05 
## ... Similar to previous best
## Run 7 stress 0.1645868 
## ... Procrustes: rmse 5.445529e-06  max resid 3.531163e-05 
## ... Similar to previous best
## Run 8 stress 0.1645868 
## ... Procrustes: rmse 2.938561e-06  max resid 1.758468e-05 
## ... Similar to previous best
## Run 9 stress 0.1645868 
## ... Procrustes: rmse 5.152021e-06  max resid 2.938458e-05 
## ... Similar to previous best
## Run 10 stress 0.1645868 
## ... Procrustes: rmse 3.21662e-06  max resid 1.860854e-05 
## ... Similar to previous best
## Run 11 stress 0.2444954 
## Run 12 stress 0.1645868 
## ... Procrustes: rmse 4.164196e-06  max resid 2.150385e-05 
## ... Similar to previous best
## Run 13 stress 0.1645868 
## ... New best solution
## ... Procrustes: rmse 1.422965e-06  max resid 4.641084e-06 
## ... Similar to previous best
## Run 14 stress 0.1645868 
## ... Procrustes: rmse 3.82778e-06  max resid 1.75597e-05 
## ... Similar to previous best
## Run 15 stress 0.2249074 
## Run 16 stress 0.1645868 
## ... Procrustes: rmse 3.358032e-06  max resid 2.109151e-05 
## ... Similar to previous best
## Run 17 stress 0.1645868 
## ... Procrustes: rmse 5.317784e-06  max resid 2.250939e-05 
## ... Similar to previous best
## Run 18 stress 0.1645868 
## ... Procrustes: rmse 1.990947e-06  max resid 7.963673e-06 
## ... Similar to previous best
## Run 19 stress 0.1645868 
## ... Procrustes: rmse 7.68143e-06  max resid 4.757098e-05 
## ... Similar to previous best
## Run 20 stress 0.1645868 
## ... Procrustes: rmse 2.829496e-06  max resid 1.629448e-05 
## ... Similar to previous best
## *** Best solution repeated 7 times
# Plot the PCoA
wUNI_station_nmds_combined_season <- 
  plot_ordination(
    physeq = scaled_rooted_physeq,
    ordination = scaled_wUNI_nmds_combined,
    color = "lakesite",
    shape = "season",
    title = "Weighted Unifrac NMDS combined season") +
  geom_point(size=5, alpha = 0.5, aes(color = lakesite)) +
  scale_shape_manual(values = c(15,16,17)) +
  scale_color_manual(values = c(lakesite_colors)) + 
  theme_bw()

wUNI_station_nmds_combined_limnion <- 
  plot_ordination(
    physeq = scaled_rooted_physeq,
    ordination = scaled_wUNI_nmds_combined,
    color = "lakesite",
    shape = "limnion",
    title = "Weighted Unifrac NMDS combined limnion") +
  geom_point(size=5, alpha = 0.5, aes(color = lakesite)) +
  scale_shape_manual(values = c(15,16,17)) +
  scale_color_manual(values = c(lakesite_colors)) + 
  theme_bw()

wUNI_station_nmds_DNA_season <- 
  plot_ordination(
    physeq = scaled_rooted_physeq_DNA,
    ordination = scaled_wUNI_nmds_DNA,
    color = "lakesite",
    shape = "season",
    title = "Weighted Unifrac NMDS DNA season") +
  geom_point(size=5, alpha = 0.5, aes(color = lakesite)) +
  scale_shape_manual(values = c(15,16,17)) +
  scale_color_manual(values = c(lakesite_colors)) + 
  theme_bw()

wUNI_station_nmds_DNA_limnion <- 
  plot_ordination(
    physeq = scaled_rooted_physeq_DNA,
    ordination = scaled_wUNI_nmds_DNA,
    color = "lakesite",
    shape = "limnion",
    title = "Weighted Unifrac NMDS DNA limnion") +
  geom_point(size=5, alpha = 0.5, aes(color = lakesite)) +
  scale_shape_manual(values = c(15,16,17)) +
  scale_color_manual(values = c(lakesite_colors)) + 
  theme_bw()

wUNI_station_nmds_RNA_season <- 
  plot_ordination(
    physeq = scaled_rooted_physeq_RNA,
    ordination = scaled_wUNI_nmds_RNA,
    color = "lakesite",
    shape = "season",
    title = "Weighted Unifrac NMDS RNA season") +
  geom_point(size=5, alpha = 0.5, aes(color = lakesite)) +
  scale_shape_manual(values = c(15,16,17)) +
  scale_color_manual(values = c(lakesite_colors)) + 
  theme_bw()

wUNI_station_nmds_RNA_limnion <- 
  plot_ordination(
    physeq = scaled_rooted_physeq_RNA,
    ordination = scaled_wUNI_nmds_RNA,
    color = "lakesite",
    shape = "limnion",
    title = "Weighted Unifrac NMDS RNA limnion") +
  geom_point(size=5, alpha = 0.5, aes(color = lakesite)) +
  scale_shape_manual(values = c(15,16,17)) +
  scale_color_manual(values = c(lakesite_colors)) + 
  theme_bw()


grid.arrange(wUNI_station_nmds_combined_season, wUNI_station_nmds_combined_limnion, ncol=2)

grid.arrange(wUNI_station_nmds_DNA_season, wUNI_station_nmds_RNA_season,ncol=2)

grid.arrange(wUNI_station_nmds_DNA_limnion, wUNI_station_nmds_RNA_limnion,ncol=2)

These NMDS plots don’t reveal any salient features, other than clustering based on lakesite.

##Statistical Significance Testing

# Calculate all three of the distance matrices

scaled_sorensen_dist_combined <- phyloseq::distance(scaled_rooted_physeq, method = "bray", binary = TRUE)

scaled_sorensen_dist_DNA <- phyloseq::distance(scaled_rooted_physeq_DNA, method = "bray", binary = TRUE)

scaled_sorensen_dist_RNA <- phyloseq::distance(scaled_rooted_physeq_RNA, method = "bray", binary = TRUE)

scaled_bray_dist_combined <- phyloseq::distance(scaled_rooted_physeq, method = "bray")

scaled_bray_dist_DNA <- phyloseq::distance(scaled_rooted_physeq_DNA, method = "bray")

scaled_bray_dist_RNA <- phyloseq::distance(scaled_rooted_physeq_RNA, method = "bray")

scaled_wUnifrac_dist_combined <- phyloseq::distance(scaled_rooted_physeq, method = "wunifrac")

scaled_wUnifrac_dist_DNA <- phyloseq::distance(scaled_rooted_physeq_DNA, method = "wunifrac")

scaled_wUnifrac_dist_RNA <- phyloseq::distance(scaled_rooted_physeq_RNA, method = "wunifrac")

# make a data frame from the sample_data
# All distance matrices will be the same metadata because they 
# originate from the same phyloseq object. 
metadata_combined <- data.frame(sample_data(scaled_rooted_physeq))
metadata_DNA <- data.frame(sample_data(scaled_rooted_physeq_DNA))
metadata_RNA <- data.frame(sample_data(scaled_rooted_physeq_RNA))

# Adonis test
# In this example we are testing the hypothesis that the five stations
# that were collected have different centroids in the ordination space 
# for each of the dissimilarity metrics, we are using a discrete variable 
adonis2(scaled_sorensen_dist_combined ~ lakesite, data = metadata_combined)
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = scaled_sorensen_dist_combined ~ lakesite, data = metadata_combined)
##           Df SumOfSqs      R2      F Pr(>F)    
## lakesite   3   5.8697 0.23886 10.983  0.001 ***
## Residual 105  18.7043 0.76114                  
## Total    108  24.5740 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis2(scaled_bray_dist_combined ~ lakesite, data = metadata_combined)
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = scaled_bray_dist_combined ~ lakesite, data = metadata_combined)
##           Df SumOfSqs      R2      F Pr(>F)    
## lakesite   3   4.4486 0.16638 6.9855  0.001 ***
## Residual 105  22.2892 0.83362                  
## Total    108  26.7379 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis2(scaled_wUnifrac_dist_combined ~ lakesite, data = metadata_combined)
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = scaled_wUnifrac_dist_combined ~ lakesite, data = metadata_combined)
##           Df SumOfSqs      R2      F Pr(>F)   
## lakesite   3   0.5782 0.08178 3.1171  0.002 **
## Residual 105   6.4923 0.91822                 
## Total    108   7.0705 1.00000                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis2(scaled_sorensen_dist_DNA ~ lakesite, data = metadata_DNA)
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = scaled_sorensen_dist_DNA ~ lakesite, data = metadata_DNA)
##          Df SumOfSqs      R2      F Pr(>F)    
## lakesite  3   3.2802 0.32063 8.1804  0.001 ***
## Residual 52   6.9503 0.67937                  
## Total    55  10.2305 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis2(scaled_bray_dist_DNA ~ lakesite, data = metadata_DNA)
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = scaled_bray_dist_DNA ~ lakesite, data = metadata_DNA)
##          Df SumOfSqs      R2      F Pr(>F)    
## lakesite  3   2.5944 0.30663 7.6653  0.001 ***
## Residual 52   5.8667 0.69337                  
## Total    55   8.4611 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis2(scaled_wUnifrac_dist_DNA ~ lakesite, data = metadata_DNA)
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = scaled_wUnifrac_dist_DNA ~ lakesite, data = metadata_DNA)
##          Df SumOfSqs      R2      F Pr(>F)    
## lakesite  3  0.28642 0.17106 3.5768  0.001 ***
## Residual 52  1.38801 0.82894                  
## Total    55  1.67444 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis2(scaled_sorensen_dist_RNA ~ lakesite, data = metadata_RNA)
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = scaled_sorensen_dist_RNA ~ lakesite, data = metadata_RNA)
##          Df SumOfSqs      R2      F Pr(>F)    
## lakesite  3   3.1514 0.25969 5.7294  0.001 ***
## Residual 49   8.9840 0.74031                  
## Total    52  12.1354 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis2(scaled_bray_dist_RNA ~ lakesite, data = metadata_RNA)
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = scaled_bray_dist_RNA ~ lakesite, data = metadata_RNA)
##          Df SumOfSqs      R2      F Pr(>F)    
## lakesite  3   2.7364 0.23186 4.9302  0.001 ***
## Residual 49   9.0656 0.76814                  
## Total    52  11.8020 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis2(scaled_wUnifrac_dist_RNA ~ lakesite, data = metadata_RNA)
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = scaled_wUnifrac_dist_RNA ~ lakesite, data = metadata_RNA)
##          Df SumOfSqs      R2      F Pr(>F)    
## lakesite  3    0.420 0.18049 3.5973  0.001 ***
## Residual 49    1.907 0.81951                  
## Total    52    2.327 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Using lakesite as our primary explanatory variable, we see that we see statistical significance for all our dissimilarity metrics when we conduct a PERMANOVA analysis. This indicates that the community composition is different at a statistically significant level between each of the different depths that were sampled and this is true for the total taxa (indicated by the DNA samples) as well as the active taxa (indicated by the RNA samples).

##BetaDispR

# Homogeneity of Disperson test with beta dispr
# Sorensen 
beta_soren_station_combined <- betadisper(scaled_sorensen_dist_combined, metadata_combined$lakesite)
permutest(beta_soren_station_combined)
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##            Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)   
## Groups      3 0.14349 0.047831 5.6975    999  0.004 **
## Residuals 105 0.88148 0.008395                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Bray-curtis 
beta_bray_station_combined <- betadisper(scaled_bray_dist_combined, metadata_combined$lakesite)
permutest(beta_bray_station_combined)
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##            Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
## Groups      3 0.04498 0.014993 1.7103    999  0.153
## Residuals 105 0.92043 0.008766
# Weighted Unifrac 
beta_unifrac_station_combined <- betadisper(scaled_wUnifrac_dist_combined, metadata_combined$lakesite)
permutest(beta_unifrac_station_combined)
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##            Df   Sum Sq    Mean Sq      F N.Perm Pr(>F)
## Groups      3 0.001152 0.00038413 0.1343    999  0.928
## Residuals 105 0.300301 0.00286001
#Homogeneity of Dispersion test with beta dispr 
#Soerson DNA
beta_soren_station_DNA <- betadisper(scaled_sorensen_dist_DNA, metadata_DNA$lakesite)
permutest(beta_soren_station_DNA)
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)  
## Groups     3 0.06805 0.0226834 2.5481    999  0.056 .
## Residuals 52 0.46291 0.0089021                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Bray-curtis DNA
beta_bray_station_DNA <- betadisper(scaled_bray_dist_DNA, metadata_DNA$lakesite)
permutest(beta_bray_station_DNA)
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
## Groups     3 0.03891 0.0129706 1.6156    999  0.217
## Residuals 52 0.41748 0.0080285
# Weighted Unifrac DNA
beta_bray_unifrac_DNA <- betadisper(scaled_wUnifrac_dist_DNA, metadata_DNA$lakesite)
permutest(beta_bray_unifrac_DNA)
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
## Groups     3 0.007897 0.0026323 1.1012    999  0.364
## Residuals 52 0.124307 0.0023905
#Homogeneity of Dispersion test with beta dispr 
#Soerson RNA
beta_soren_station_RNA <- betadisper(scaled_sorensen_dist_RNA, metadata_RNA$lakesite)
permutest(beta_soren_station_RNA)
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
## Groups     3 0.09634 0.032114 3.6571    999   0.02 *
## Residuals 49 0.43028 0.008781                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Bray-curtis RNA
beta_bray_station_RNA <- betadisper(scaled_bray_dist_RNA, metadata_RNA$lakesite)
permutest(beta_bray_station_RNA)
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
## Groups     3 0.07377 0.024588 2.3909    999  0.096 .
## Residuals 49 0.50393 0.010284                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Weighted Unifrac rNA
beta_bray_unifrac_RNA <- betadisper(scaled_wUnifrac_dist_RNA, metadata_RNA$lakesite)
permutest(beta_bray_unifrac_RNA)
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df   Sum Sq    Mean Sq      F N.Perm Pr(>F)
## Groups     3 0.002733 0.00091112 0.2955    999  0.835
## Residuals 49 0.151078 0.00308323

The high p-values from the betadisper test indicate that the variability is fairly homogeneous across our different samples.

Taxonomic Composition

#Phylum

# Set the phylum colors
phylum_colors <- c(
  Acidobacteriota = "navy", 
  Actinobacteriota = "darkslategray2", 
  Armatimonadota = "deeppink1",
  Alphaproteobacteria = "plum2", 
  Bacteroidota = "gold", 
  Betaproteobacteria = "plum1", 
  Bdellovibrionota = "red1",
  Chloroflexi="black", 
  Crenarchaeota = "firebrick",
  Cyanobacteria = "limegreen",
  Deltaproteobacteria = "grey", 
  Desulfobacterota="magenta",
  Firmicutes = "#3E9B96",
  Gammaproteobacteria = "greenyellow",
  "Marinimicrobia (SAR406 clade)" = "yellow",
  Myxococcota = "#B5D6AA",
  Nitrospirota = "palevioletred1",
  Proteobacteria = "royalblue",
  Planctomycetota = "darkorange", 
  "SAR324 clade(Marine group B)" = "olivedrab",
  #Proteobacteria_unclassified = "greenyellow",
  Thermoplasmatota = "green",
  Verrucomicrobiota = "darkorchid1")
 # Other = "grey")
# Calculate the phylum relative abundance 
# Make different dfs for combined, DNA and RNA samples
phylum_df_combined <- 
  scaled_rooted_physeq %>%
  # agglomerate at the phylum level 
  tax_glom(taxrank = "Phylum") %>% 
  # Transform counts to relative abundance 
  transform_sample_counts(function (x) {x/sum(x)}) %>%
  # Melt to a long format 
  psmelt() %>%
  # Filter out Phyla that are <1% - get rid of low abundant Phyla
  dplyr::filter(Abundance > 0.01)

phylum_df_DNA <- 
  scaled_rooted_physeq_DNA %>%
  # agglomerate at the phylum level 
  tax_glom(taxrank = "Phylum") %>% 
  # Transform counts to relative abundance 
  transform_sample_counts(function (x) {x/sum(x)}) %>%
  # Melt to a long format 
  psmelt() %>%
  # Filter out Phyla that are <1% - get rid of low abundant Phyla
  dplyr::filter(Abundance > 0.01)

phylum_df_RNA <- 
  scaled_rooted_physeq_RNA %>%
  # agglomerate at the phylum level 
  tax_glom(taxrank = "Phylum") %>% 
  # Transform counts to relative abundance 
  transform_sample_counts(function (x) {x/sum(x)}) %>%
  # Melt to a long format 
  psmelt() %>%
  # Filter out Phyla that are <1% - get rid of low abundant Phyla
  dplyr::filter(Abundance > 0.01)

#Combined Phyla Plots

# Stacked Bar Plot With All phyla 
# Plot Phylum Abundances - make sure to load phylum_colors 
phylum_df_combined %>%
  # Warning: It's important to have one sample per x value, 
  # otherwise, it will take the sum between multiple samples
  dplyr::filter(limnion == "Bottom") %>% 
  ggplot(aes(x = lakesite, y = Abundance, fill = Phylum)) + 
  facet_grid(Phylum~season, scale = "free") +
  geom_bar(stat = "identity", color = "black") + 
  labs(title = "Bottom Limnion Phylum Composition") + 
  scale_fill_manual(values = phylum_colors) + 
  theme_bw() + 
  theme(axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1))

phylum_df_combined %>%
  # Warning: It's important to have one sample per x value, 
  # otherwise, it will take the sum between multiple samples
  dplyr::filter(limnion == "Middle") %>% 
  ggplot(aes(x = lakesite, y = Abundance, fill = Phylum)) + 
  facet_grid(Phylum~season, scale = "free") +
  geom_bar(stat = "identity", color = "black") + 
  labs(title = "Middle Limnion Phylum Composition") + 
  scale_fill_manual(values = phylum_colors) + 
  theme_bw() + 
  theme(axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1))

phylum_df_combined %>%
  # Warning: It's important to have one sample per x value, 
  # otherwise, it will take the sum between multiple samples
  dplyr::filter(limnion == "Top") %>% 
  ggplot(aes(x = lakesite, y = Abundance, fill = Phylum)) + 
  facet_grid(Phylum~season, scale = "free") +
  geom_bar(stat = "identity", color = "black") + 
  labs(title = "Top Limnion Phylum Composition") + 
  scale_fill_manual(values = phylum_colors) + 
  theme_bw() + 
  theme(axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1))

#Phyla plots from DNA samples

phylum_df_DNA %>%
  # Warning: It's important to have one sample per x value, 
  # otherwise, it will take the sum between multiple samples
  dplyr::filter(limnion == "Top") %>% 
  ggplot(aes(x = lakesite, y = Abundance, fill = Phylum)) + 
  facet_grid(Phylum~season, scale = "free") +
  geom_bar(stat = "identity", color = "black") + 
  labs(title = "Top Limnion Phylum Composition DNA") + 
  scale_fill_manual(values = phylum_colors) + 
  theme_bw() + 
  theme(axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1))

phylum_df_DNA %>%
  # Warning: It's important to have one sample per x value, 
  # otherwise, it will take the sum between multiple samples
  dplyr::filter(limnion == "Middle") %>% 
  ggplot(aes(x = lakesite, y = Abundance, fill = Phylum)) + 
  facet_grid(Phylum~season, scale = "free") +
  geom_bar(stat = "identity", color = "black") + 
  labs(title = "Middle Limnion Phylum Composition DNA") + 
  scale_fill_manual(values = phylum_colors) + 
  theme_bw() + 
  theme(axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1))

phylum_df_DNA %>%
  # Warning: It's important to have one sample per x value, 
  # otherwise, it will take the sum between multiple samples
  dplyr::filter(limnion == "Bottom") %>% 
  ggplot(aes(x = lakesite, y = Abundance, fill = Phylum)) + 
  facet_grid(Phylum~season, scale = "free") +
  geom_bar(stat = "identity", color = "black") + 
  labs(title = "Bottom Limnion Phylum Composition DNA") + 
  scale_fill_manual(values = phylum_colors) + 
  theme_bw() + 
  theme(axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1))

Top Limnion: For the total taxa based on lakesite we see slight increase in Actinobacteriota in summer throughout each lakesite and this seems to be the most abundant phyla. These are not the dominant active taxa, however. The dominant taxa according to the RNA samples are Verrucomicrobiota and Proteobacteria.Armatiomondota seem primarily to be at the surface and do not make appear in the active taxa. Cyanobacteria seem to be quite abundant in fall and summer across lake site. Interestingly the 15 meter lakesite seems to be unique in having cyanobateria in spring.

Middle Limnion: The only middle limnion lakesite sampled was 110m from the shore. The total taxa for the middle limnion is dominated by actinobacteriota and bacteroidota across season with marked increase in verrucomicrobiota in summer. This tracks well with the RNA samples, indicating similar trends between the total and active taxa.

Bottom Limnion: We similar trends of actinobacteriota to the top and middle limnion with similar trends across season. A point of note is that there seems to be no cyanobacteria present at any lakesite in spring. This limnion is also dominated by verrucomicriobiota. The active taxa track nicely with the total taxa for this depth, with the exception of Chloroflexi that make up a sizeable portion of the active taxa.

#Phyla from RNA samples

phylum_df_RNA <- 
  scaled_rooted_physeq_RNA %>%
  # agglomerate at the phylum level 
  tax_glom(taxrank = "Phylum") %>% 
  # Transform counts to relative abundance 
  transform_sample_counts(function (x) {x/sum(x)}) %>%
  # Melt to a long format 
  psmelt() %>%
  # Filter out Phyla that are <1% - get rid of low abundant Phyla
  dplyr::filter(Abundance > 0.01)

phylum_df_RNA %>%
  # Warning: It's important to have one sample per x value, 
  # otherwise, it will take the sum between multiple samples
  dplyr::filter(limnion == "Middle") %>% 
  ggplot(aes(x = lakesite, y = Abundance, fill = Phylum)) + 
  facet_grid(Phylum~season, scale = "free") +
  geom_bar(stat = "identity", color = "black") + 
  labs(title = "Middle Limnion Phylum Composition RNA") + 
  scale_fill_manual(values = phylum_colors) + 
  theme_bw() + 
  theme(axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1)) 

phylum_df_RNA %>%
  # Warning: It's important to have one sample per x value, 
  # otherwise, it will take the sum between multiple samples
  dplyr::filter(limnion == "Bottom") %>% 
  ggplot(aes(x = lakesite, y = Abundance, fill = Phylum)) + 
  facet_grid(Phylum~season, scale = "free") +
  geom_bar(stat = "identity", color = "black") + 
  labs(title = "Bottom Limnion Phylum Composition RNA") + 
  scale_fill_manual(values = phylum_colors) + 
  theme_bw() + 
  theme(axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1))

phylum_df_RNA %>%
  # Warning: It's important to have one sample per x value, 
  # otherwise, it will take the sum between multiple samples
  dplyr::filter(limnion == "Top") %>% 
  ggplot(aes(x = lakesite, y = Abundance, fill = Phylum)) + 
  facet_grid(Phylum~season, scale = "free") +
  geom_bar(stat = "identity", color = "black") + 
  labs(title = "Top Limnion Phylum Composition RNA") + 
  scale_fill_manual(values = phylum_colors) + 
  theme_bw() + 
  theme(axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1))

##Family

# Calculate the Family relative abundance for each sample
# # Make separate dfs for combined, DNA and RNA samples
family_df_combined <- 
  scaled_rooted_physeq %>%
  # agglomerate at the phylum level 
  tax_glom(taxrank = "Family") %>% 
  # Transform counts to relative abundance 
  transform_sample_counts(function (x) {x/sum(x)}) %>%
  # Melt to a long format 
  psmelt() %>%
  # Filter out Phyla that are <1% - get rid of low abundant Phyla
  dplyr::filter(Abundance > 0.01) 

family_df_DNA <- 
  scaled_rooted_physeq_DNA %>%
  # agglomerate at the phylum level 
  tax_glom(taxrank = "Family") %>% 
  # Transform counts to relative abundance 
  transform_sample_counts(function (x) {x/sum(x)}) %>%
  # Melt to a long format 
  psmelt() %>%
  # Filter out Phyla that are <1% - get rid of low abundant Phyla
  dplyr::filter(Abundance > 0.01) 

family_df_RNA <- 
  scaled_rooted_physeq_RNA %>%
  # agglomerate at the phylum level 
  tax_glom(taxrank = "Family") %>% 
  # Transform counts to relative abundance 
  transform_sample_counts(function (x) {x/sum(x)}) %>%
  # Melt to a long format 
  psmelt() %>%
  # Filter out Phyla that are <1% - get rid of low abundant Phyla
  dplyr::filter(Abundance > 0.01) 
#Judging from the Phylum plots we identified Cyanobacteria,Actinobacteriota and Verrucomicrobiota as the most interesting Phyla to discet further as they seem to be quite abudant and vary between season and limnion.
# Cyanobacteria Families in DNA sample 
# Plot family 
cyano_season_DNA <- family_df_DNA %>%
  dplyr::filter(Phylum == "Cyanobacteria") %>%
  # build the plot 
  ggplot(aes(x = lakesite, y = Abundance, 
             fill = lakesite, color = lakesite)) + 
  facet_grid(Family~season, scale = "free") +
  geom_boxplot(alpha = 0.5, outlier.shape = NA) + # outliers not plotted here in boxplot 
  geom_jitter() + 
  theme_bw() + 
  labs(title = "Cyanobacteria Family Abundance") + 
  scale_color_manual(values = lakesite_colors) + 
  scale_fill_manual(values = lakesite_colors) + 
  theme(axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1),
        legend.position = "bottom")

actin_season_dna <- family_df_DNA %>%
  dplyr::filter(Phylum == "Actinobacteriota") %>%
  # build the plot 
  ggplot(aes(x = lakesite, y = Abundance, 
             fill = lakesite, color = lakesite)) + 
  facet_grid(Family~season, scale = "free") +
  geom_boxplot(alpha = 0.5, outlier.shape = NA) + # outliers not plotted here in boxplot 
  geom_jitter() + 
  theme_bw() + 
  labs(title = "Actinobacteriota Family Abundance") + 
  scale_color_manual(values = lakesite_colors) + 
  scale_fill_manual(values = lakesite_colors) + 
  theme(axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1),
        legend.position = "bottom")
#actin_dna sanity check

verru_season_dna <- family_df_DNA %>%
  dplyr::filter(Phylum == "Verrucomicrobiota") %>%
  # build the plot 
  ggplot(aes(x = lakesite, y = Abundance, 
             fill = lakesite, color = lakesite)) + 
  facet_grid(Family~season, scale = "free") +
  geom_boxplot(alpha = 0.5, outlier.shape = NA) + # outliers not plotted here in boxplot 
  geom_jitter() + 
  theme_bw() + 
  labs(title = " VerrucomicrobiotaFamily Abundance") + 
  scale_color_manual(values = lakesite_colors) + 
  scale_fill_manual(values = lakesite_colors) + 
  theme(axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1),
        legend.position = "bottom")
#verru_dna sanity check

#Now we will check the DNA families faceting by limnion
# Cyanobacteria Families in DNA sample 
# Plot family 
cyano_lim_DNA <- family_df_DNA %>%
  dplyr::filter(Phylum == "Cyanobacteria") %>%
  # build the plot 
  ggplot(aes(x = lakesite, y = Abundance, 
             fill = lakesite, color = lakesite)) + 
  facet_grid(Family~limnion, scale = "free") +
  geom_boxplot(alpha = 0.5, outlier.shape = NA) + # outliers not plotted here in boxplot 
  geom_jitter() + 
  theme_bw() + 
  labs(title = "Cyanobacteria Family Abundance") + 
  scale_color_manual(values = lakesite_colors) + 
  scale_fill_manual(values = lakesite_colors) + 
  theme(axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1),
        legend.position = "bottom")

actin_lim_dna <- family_df_DNA %>%
  dplyr::filter(Phylum == "Actinobacteriota") %>%
  # build the plot 
  ggplot(aes(x = lakesite, y = Abundance, 
             fill = lakesite, color = lakesite)) + 
  facet_grid(Family~limnion, scale = "free") +
  geom_boxplot(alpha = 0.5, outlier.shape = NA) + # outliers not plotted here in boxplot 
  geom_jitter() + 
  theme_bw() + 
  labs(title = "Actinobacteriota Family Abundance") + 
  scale_color_manual(values = lakesite_colors) + 
  scale_fill_manual(values = lakesite_colors) + 
  theme(axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1),
        legend.position = "bottom")
#actin_dna sanity check

verru_limnion_dna <- family_df_DNA %>%
  dplyr::filter(Phylum == "Verrucomicrobiota") %>%
  # build the plot 
  ggplot(aes(x = lakesite, y = Abundance, 
             fill = lakesite, color = lakesite)) + 
  facet_grid(Family~limnion, scale = "free") +
  geom_boxplot(alpha = 0.5, outlier.shape = NA) + # outliers not plotted here in boxplot 
  geom_jitter() + 
  theme_bw() + 
  labs(title = " VerrucomicrobiotaFamily Abundance") + 
  scale_color_manual(values = lakesite_colors) + 
  scale_fill_manual(values = lakesite_colors) + 
  theme(axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1),
        legend.position = "bottom")
#verru_dna sanity check

#RNA Samples Plot

cyano_season_RNA <- family_df_RNA %>%
  dplyr::filter(Phylum == "Cyanobacteria") %>%
  # build the plot 
  ggplot(aes(x = lakesite, y = Abundance, 
             fill = lakesite, color = lakesite)) + 
  facet_grid(Family~season, scale = "free") +
  geom_boxplot(alpha = 0.5, outlier.shape = NA) + # outliers not plotted here in boxplot 
  geom_jitter() + 
  theme_bw() + 
  labs(title = "Cyanobacteria RNA Family Abundance") + 
  scale_color_manual(values = lakesite_colors) + 
  scale_fill_manual(values = lakesite_colors) + 
  theme(axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1),
        legend.position = "bottom")

actin_season_rna <- family_df_RNA %>%
  dplyr::filter(Phylum == "Actinobacteriota") %>%
  # build the plot 
  ggplot(aes(x = lakesite, y = Abundance, 
             fill = lakesite, color = lakesite)) + 
  facet_grid(Family~season, scale = "free") +
  geom_boxplot(alpha = 0.5, outlier.shape = NA) + # outliers not plotted here in boxplot 
  geom_jitter() + 
  theme_bw() + 
  labs(title = "Actinobacteriota Family RNA Abundance") + 
  scale_color_manual(values = lakesite_colors) + 
  scale_fill_manual(values = lakesite_colors) + 
  theme(axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1),
        legend.position = "bottom")
#actin_dna sanity check

verru_season_rna <- family_df_RNA %>%
  dplyr::filter(Phylum == "Verrucomicrobiota") %>%
  # build the plot 
  ggplot(aes(x = lakesite, y = Abundance, 
             fill = lakesite, color = lakesite)) + 
  facet_grid(Family~season, scale = "free") +
  geom_boxplot(alpha = 0.5, outlier.shape = NA) + # outliers not plotted here in boxplot 
  geom_jitter() + 
  theme_bw() + 
  labs(title = " Verrucomicrobiota Family RNA Abundance") + 
  scale_color_manual(values = lakesite_colors) + 
  scale_fill_manual(values = lakesite_colors) + 
  theme(axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1),
        legend.position = "bottom")
#RNA Limnion Samples
cyano_lim_RNA <- family_df_RNA %>%
  dplyr::filter(Phylum == "Cyanobacteria") %>%
  # build the plot 
  ggplot(aes(x = lakesite, y = Abundance, 
             fill = lakesite, color = lakesite)) + 
  facet_grid(Family~limnion, scale = "free") +
  geom_boxplot(alpha = 0.5, outlier.shape = NA) + # outliers not plotted here in boxplot 
  geom_jitter() + 
  theme_bw() + 
  labs(title = "Cyanobacteria RNA Family Abundance") + 
  scale_color_manual(values = lakesite_colors) + 
  scale_fill_manual(values = lakesite_colors) + 
  theme(axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1),
        legend.position = "bottom")

actin_lim_rna <- family_df_RNA %>%
  dplyr::filter(Phylum == "Actinobacteriota") %>%
  # build the plot 
  ggplot(aes(x = lakesite, y = Abundance, 
             fill = lakesite, color = lakesite)) + 
  facet_grid(Family~limnion, scale = "free") +
  geom_boxplot(alpha = 0.5, outlier.shape = NA) + # outliers not plotted here in boxplot 
  geom_jitter() + 
  theme_bw() + 
  labs(title = "Actinobacteriota Family RNA Abundance") + 
  scale_color_manual(values = lakesite_colors) + 
  scale_fill_manual(values = lakesite_colors) + 
  theme(axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1),
        legend.position = "bottom")
#actin_dna sanity check

verru_lim_rna <- family_df_RNA %>%
  dplyr::filter(Phylum == "Verrucomicrobiota") %>%
  # build the plot 
  ggplot(aes(x = lakesite, y = Abundance, 
             fill = lakesite, color = lakesite)) + 
  facet_grid(Family~limnion, scale = "free") +
  geom_boxplot(alpha = 0.5, outlier.shape = NA) + # outliers not plotted here in boxplot 
  geom_jitter() + 
  theme_bw() + 
  labs(title = " Verrucomicrobiota Family RNA Abundance") + 
  scale_color_manual(values = lakesite_colors) + 
  scale_fill_manual(values = lakesite_colors) + 
  theme(axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1),
        legend.position = "bottom")

cyano_lim_DNA

cyano_season_DNA

verru_limnion_dna

verru_season_rna

actin_season_dna

actin_season_rna

cyano_lim_RNA

verru_season_rna

actin_lim_rna

cyano_season_RNA

actin_season_rna

verru_season_rna

##Session Information

devtools::session_info()
## ─ Session info ───────────────────────────────────────────────────────────────
##  setting  value
##  version  R version 4.3.2 (2023-10-31)
##  os       Rocky Linux 9.0 (Blue Onyx)
##  system   x86_64, linux-gnu
##  ui       X11
##  language (EN)
##  collate  en_US.UTF-8
##  ctype    en_US.UTF-8
##  tz       America/New_York
##  date     2024-04-30
##  pandoc   3.1.1 @ /usr/lib/rstudio-server/bin/quarto/bin/tools/ (via rmarkdown)
## 
## ─ Packages ───────────────────────────────────────────────────────────────────
##  package          * version    date (UTC) lib source
##  ade4               1.7-22     2023-02-06 [1] CRAN (R 4.3.2)
##  ape                5.7-1      2023-03-13 [2] CRAN (R 4.3.2)
##  Biobase            2.62.0     2023-10-24 [2] Bioconductor
##  BiocGenerics       0.48.1     2023-11-01 [2] Bioconductor
##  biomformat         1.30.0     2023-10-24 [1] Bioconductor
##  Biostrings         2.70.1     2023-10-25 [2] Bioconductor
##  bitops             1.0-7      2021-04-24 [2] CRAN (R 4.3.2)
##  bslib              0.5.1      2023-08-11 [2] CRAN (R 4.3.2)
##  cachem             1.0.8      2023-05-01 [2] CRAN (R 4.3.2)
##  callr              3.7.3      2022-11-02 [2] CRAN (R 4.3.2)
##  cli                3.6.1      2023-03-23 [2] CRAN (R 4.3.2)
##  cluster            2.1.4      2022-08-22 [2] CRAN (R 4.3.2)
##  codetools          0.2-19     2023-02-01 [2] CRAN (R 4.3.2)
##  colorspace         2.1-0      2023-01-23 [2] CRAN (R 4.3.2)
##  crayon             1.5.2      2022-09-29 [2] CRAN (R 4.3.2)
##  data.table         1.14.8     2023-02-17 [2] CRAN (R 4.3.2)
##  devtools         * 2.4.4      2022-07-20 [2] CRAN (R 4.2.1)
##  digest             0.6.33     2023-07-07 [2] CRAN (R 4.3.2)
##  dplyr            * 1.1.3      2023-09-03 [2] CRAN (R 4.3.2)
##  ellipsis           0.3.2      2021-04-29 [2] CRAN (R 4.3.2)
##  evaluate           0.23       2023-11-01 [2] CRAN (R 4.3.2)
##  fansi              1.0.5      2023-10-08 [2] CRAN (R 4.3.2)
##  farver             2.1.1      2022-07-06 [2] CRAN (R 4.3.2)
##  fastmap            1.1.1      2023-02-24 [2] CRAN (R 4.3.2)
##  forcats          * 1.0.0      2023-01-29 [1] CRAN (R 4.3.2)
##  foreach            1.5.2      2022-02-02 [2] CRAN (R 4.3.2)
##  fs                 1.6.3      2023-07-20 [2] CRAN (R 4.3.2)
##  generics           0.1.3      2022-07-05 [2] CRAN (R 4.3.2)
##  GenomeInfoDb       1.38.0     2023-10-24 [2] Bioconductor
##  GenomeInfoDbData   1.2.11     2023-11-07 [2] Bioconductor
##  ggplot2          * 3.5.0      2024-02-23 [2] CRAN (R 4.3.2)
##  glue               1.6.2      2022-02-24 [2] CRAN (R 4.3.2)
##  gridExtra        * 2.3        2017-09-09 [1] CRAN (R 4.3.2)
##  gtable             0.3.4      2023-08-21 [2] CRAN (R 4.3.2)
##  highr              0.10       2022-12-22 [2] CRAN (R 4.3.2)
##  hms                1.1.3      2023-03-21 [1] CRAN (R 4.3.2)
##  htmltools          0.5.7      2023-11-03 [2] CRAN (R 4.3.2)
##  htmlwidgets        1.6.2      2023-03-17 [2] CRAN (R 4.3.2)
##  httpuv             1.6.12     2023-10-23 [2] CRAN (R 4.3.2)
##  igraph             1.5.1      2023-08-10 [2] CRAN (R 4.3.2)
##  IRanges            2.36.0     2023-10-24 [2] Bioconductor
##  iterators          1.0.14     2022-02-05 [2] CRAN (R 4.3.2)
##  jquerylib          0.1.4      2021-04-26 [2] CRAN (R 4.3.2)
##  jsonlite           1.8.7      2023-06-29 [2] CRAN (R 4.3.2)
##  knitr              1.45       2023-10-30 [2] CRAN (R 4.3.2)
##  labeling           0.4.3      2023-08-29 [2] CRAN (R 4.3.2)
##  later              1.3.1      2023-05-02 [2] CRAN (R 4.3.2)
##  lattice          * 0.21-9     2023-10-01 [2] CRAN (R 4.3.2)
##  lifecycle          1.0.3      2022-10-07 [2] CRAN (R 4.3.2)
##  lubridate        * 1.9.3      2023-09-27 [1] CRAN (R 4.3.2)
##  magrittr           2.0.3      2022-03-30 [2] CRAN (R 4.3.2)
##  MASS               7.3-60     2023-05-04 [2] CRAN (R 4.3.2)
##  Matrix             1.6-1.1    2023-09-18 [2] CRAN (R 4.3.2)
##  memoise            2.0.1      2021-11-26 [2] CRAN (R 4.3.2)
##  mgcv               1.9-0      2023-07-11 [2] CRAN (R 4.3.2)
##  mime               0.12       2021-09-28 [2] CRAN (R 4.3.2)
##  miniUI             0.1.1.1    2018-05-18 [2] CRAN (R 4.3.2)
##  multtest           2.58.0     2023-10-24 [1] Bioconductor
##  munsell            0.5.0      2018-06-12 [2] CRAN (R 4.3.2)
##  nlme               3.1-163    2023-08-09 [2] CRAN (R 4.3.2)
##  pacman             0.5.1      2019-03-11 [1] CRAN (R 4.3.2)
##  patchwork        * 1.2.0.9000 2024-04-28 [1] Github (thomasp85/patchwork@d943757)
##  permute          * 0.9-7      2022-01-27 [1] CRAN (R 4.3.2)
##  phyloseq         * 1.46.0     2023-10-24 [1] Bioconductor
##  pillar             1.9.0      2023-03-22 [2] CRAN (R 4.3.2)
##  pkgbuild           1.4.2      2023-06-26 [2] CRAN (R 4.3.2)
##  pkgconfig          2.0.3      2019-09-22 [2] CRAN (R 4.3.2)
##  pkgload            1.3.3      2023-09-22 [2] CRAN (R 4.3.2)
##  plyr               1.8.9      2023-10-02 [2] CRAN (R 4.3.2)
##  prettyunits        1.2.0      2023-09-24 [2] CRAN (R 4.3.2)
##  processx           3.8.2      2023-06-30 [2] CRAN (R 4.3.2)
##  profvis            0.3.8      2023-05-02 [2] CRAN (R 4.3.2)
##  promises           1.2.1      2023-08-10 [2] CRAN (R 4.3.2)
##  ps                 1.7.5      2023-04-18 [2] CRAN (R 4.3.2)
##  purrr            * 1.0.2      2023-08-10 [2] CRAN (R 4.3.2)
##  R6                 2.5.1      2021-08-19 [2] CRAN (R 4.3.2)
##  Rcpp               1.0.11     2023-07-06 [2] CRAN (R 4.3.2)
##  RCurl              1.98-1.13  2023-11-02 [2] CRAN (R 4.3.2)
##  readr            * 2.1.5      2024-01-10 [1] CRAN (R 4.3.2)
##  remotes            2.4.2.1    2023-07-18 [2] CRAN (R 4.3.2)
##  reshape2           1.4.4      2020-04-09 [2] CRAN (R 4.3.2)
##  rhdf5              2.46.1     2023-11-29 [1] Bioconductor 3.18 (R 4.3.2)
##  rhdf5filters       1.14.1     2023-11-06 [1] Bioconductor
##  Rhdf5lib           1.24.2     2024-02-07 [1] Bioconductor 3.18 (R 4.3.2)
##  rlang              1.1.2      2023-11-04 [2] CRAN (R 4.3.2)
##  rmarkdown          2.25       2023-09-18 [2] CRAN (R 4.3.2)
##  rstudioapi         0.15.0     2023-07-07 [2] CRAN (R 4.3.2)
##  S4Vectors          0.40.1     2023-10-26 [2] Bioconductor
##  sass               0.4.7      2023-07-15 [2] CRAN (R 4.3.2)
##  scales             1.3.0      2023-11-28 [2] CRAN (R 4.3.2)
##  sessioninfo        1.2.2      2021-12-06 [2] CRAN (R 4.3.2)
##  shiny              1.7.5.1    2023-10-14 [2] CRAN (R 4.3.2)
##  stringi            1.7.12     2023-01-11 [2] CRAN (R 4.3.2)
##  stringr          * 1.5.0      2022-12-02 [2] CRAN (R 4.3.2)
##  survival           3.5-7      2023-08-14 [2] CRAN (R 4.3.2)
##  tibble           * 3.2.1      2023-03-20 [2] CRAN (R 4.3.2)
##  tidyr            * 1.3.0      2023-01-24 [2] CRAN (R 4.3.2)
##  tidyselect         1.2.0      2022-10-10 [2] CRAN (R 4.3.2)
##  tidyverse        * 2.0.0      2023-02-22 [1] CRAN (R 4.3.2)
##  timechange         0.3.0      2024-01-18 [1] CRAN (R 4.3.2)
##  tzdb               0.4.0      2023-05-12 [1] CRAN (R 4.3.2)
##  urlchecker         1.0.1      2021-11-30 [2] CRAN (R 4.3.2)
##  usethis          * 2.2.2      2023-07-06 [2] CRAN (R 4.3.2)
##  utf8               1.2.4      2023-10-22 [2] CRAN (R 4.3.2)
##  vctrs              0.6.4      2023-10-12 [2] CRAN (R 4.3.2)
##  vegan            * 2.6-4      2022-10-11 [1] CRAN (R 4.3.2)
##  withr              2.5.2      2023-10-30 [2] CRAN (R 4.3.2)
##  xfun               0.41       2023-11-01 [2] CRAN (R 4.3.2)
##  xtable             1.8-4      2019-04-21 [2] CRAN (R 4.3.2)
##  XVector            0.42.0     2023-10-24 [2] Bioconductor
##  yaml               2.3.7      2023-01-23 [2] CRAN (R 4.3.2)
##  zlibbioc           1.48.0     2023-10-24 [2] Bioconductor
## 
##  [1] /home/dt473/R/x86_64-pc-linux-gnu-library/4.3
##  [2] /programs/R-4.3.2/library
## 
## ──────────────────────────────────────────────────────────────────────────────